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Abstract. NGC 4449 is an active star-forming dwarf galaxy of Magellanic type. From radio observations, van 



Wocrden et al. (1975) found an extended Hl-halo around NG C 444 9 which is at least a factor of 10 larger than 



the optical diameter D25 ~ 5.6 kpc. Recently, Hunter et al. (1998) discerned details in the Hl-halo: a disc-like 
feature around the center of NGC 4449 and a lopsided arm structure. 

We combined several N-body methods in order to investigate the interaction scenario between NGC 4449 and 
DDO 125, a close companion in projected space. In a first step fast restricted N-body models are used to confine 
a region in parameter space reproducing the main observational features. In a second step a genetic algorithm is 
applied for a uniqueness test of our preferred parameter set. We show that our genetic algorithm reliably recovers 
orbital parameters, provided that the data are sufficiently accurate, i.e. all the key features are included. In the 
third step the results of the restricted N-body models are compared with self-consistent N-body simulations. In 
the case of NGC 4449, the applicability of the simple restricted N-body calculations is demonstrated. Additionally, 
it is shown that the HI gas can be modeled here by a purely stellar dynamical approach. 

In a series of simulations, we demonstrate that the observed features of the extended HI disc can be explained by 
a gravitational interaction between NGC 4449 and DDO 125. According to these calculations the closest approach 
between both galaxies happened ~ 4 — 6 • 10* yr ago at a minimum distance of ~ 25 kpc on a parabolic or slightly 
elliptic orbit. In the case of an encounter scenario, the dynamical mass of DDO 125 should not be smaller than 
10% of NGC 4449's mass. Before the encounter, the observed HI gas was arranged in a disc with a radius of 35--40 
kpc around the center of NGC 4449. It had the same orientation as the central ellipsoidal HI structure. The origin 
of this disc is still unclear, but it might have been caused by a previous interaction. 

Key words, galaxies: interactions - galaxies: kinematics and dynamics - individual: NGC 4449, DDO 125 - 
methods: N-body simulations - methods: data analysis 



1. Introduction 

During the last decades the classical picture of galaxies 
being Weltinseln (e.g. Kant 1755; v. Humboldt 185C), i.e. 
island universes, has completely changed. In the 1920s, 
galaxies were thought to be in general islands, i.e. isolated 
stellar systems, with just a few obvious exceptions like 
the M51 system. Systems not fitting the Hubble classifi- 
cation have been neglected - by definition - as peculiar. 
Interest in these systems has increased strongly since the 
new catalogues by Vorontsov-Velyaminov (1959) and Arp 
( |1966D became available. They demonstrated that even 
peculiar objects have common features like bridges, tails, 
or rings. Better observational instruments allowing for 
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deeper images and new wavelength ranges have revealed 
more complex structures. Galaxies formerly classified as 
non-interacting may have had some gravitational interac- 
tion in the past which is still reflected in e.g. warps, flares, 
or thick discs. Especially useful are HI observations which 
can cover a much larger radial extension than the opti- 
cal images. Therefore, HI is a very good tracer for tidal 
interactions. 

The theoretical understanding of interacting galaxies 
suffered for a long time from the lack of computational 
power allowing for a numerical solution of the gravita- 
tional N-bo dy pr oblem. After a remarkable treatment by 
Holmberg ( 1941 ) who built an analogue computer (con- 
sisting of light bulbs and photo cells) to determine the 
gravitational force, it took 20 years until N-body simu- 
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lations were performed on a general purpose computer 
(Pfleiderer & Siedentopf 1961): The basic idea of these 
restricted N-body simulations is the assumption that the 
potential of interacting galaxies can be adequately mod- 
eled by two particles which represent the galactic masses 
and move under their mutual gravitation, i.e. on Keplerian 
orbits. With these assumptions all the other particles are 
just test particles, and the complete N-body problem is 
reduced to N single body problems for a time-dependent 
potential. I n a r emarkable series of simulations Toomre 
& Toomre ( 1972 ) applied this technique to determine the 
parameters of some well-studied interacting systems like 
Arp 295, M 51 -h NGC 5195, or NGC 4038/39. 

New N-body techniques have been developed which in- 
creased the accessible particle number by many orders of 
magnitude . E.g. theTREE- method (Barnes & Hut |l986 ; 



Hernquist 1987 , 1990 ) in which the organization of the 
force calculation - the most time-consuming part of N- 
body calculations - adapts to clumpy mass distributions 
for a moderate computational price (~ 0{N log N)) com- 
pared to direct simulations (~ 0{N'^)). By this, Barnes 
( |1988| ) was able to simulate encounters of disc galaxies 
including all dynamical components, i.e. the disc, the 
bulge, and the halo as N-body systems. Compared to 
faster grid-based methods (e.g. Sellwood 1980) or expan- 
sion methods (e.g. Hernquist & Ostriker 1992) direct N- 



body simulations (or semi-direct methods like TREEs) 
are more flexible with respect to strongly varying geome- 
tries and scalelengths. An alternative to these techniques 
are special-purpose computers (such as th e machines of 
the GRAPE project (Sugimoto et al. |l99C| )). They imple- 
mented Newton's law of gravity (modified for gravitational 
softening) in the hardware which allows for a very fast di- 
rect determination of the gravitational forces (however, 
there still exists the N'^ bottleneck). E.g. for simulations 
with N = 10^ particles a GRAPE3af (with 8 GRAPE 
processors) is competitive with a TREE-code running on 
a CRAY T90. 

Two main methods (based on particles) to treat a dissi- 
pative, star-forming gas have been applied: The smoothed 
particle hydrodynamics scheme (SPH) solves the gasdy- 
namical equations, by this emphasizing the diffuse nature 
of the interstellar medium (ISM) (e.g. Hernquist & Katz 



198£). An alternative approach focuses on the dumpiness 



of the ISM treated as sticky particles: Without physical 
collisions or close encounters of clouds they move in bal- 
listic orbits like stars. However, in the case of a collision 
the clouds might merge or lose kinetic energy depending 



on the adopted microphysics (e.g. Casoli & Combes 1982 
Palous et al. |1993| Theis & Hensler |1993D. 



Relatively few papers exist modeling special objects. 
This deficiency is the result of two factors: First high res- 
olution data in configuration and velocity space are re- 
quired, covering a large fraction of the space between the 
interacting galaxies. In principle, HI would be suitable, 
however, there are just a few observational sites which 
give data of sufficient quality. 



The second problem is the large parameter space for 
a galactic interaction resulting in two connected difficul- 
ties: finding a good fit and determining its uniqueness (or 
other acceptable parameter sets). Observationally, only 
three kinematical quantities - the projected position on 
the sky and the line-of-sight velocity - can be measured. 
Another parameter, the galactic mass, depends on the 
availability of velocity data, the determination of the dis- 
tance, and the reliability of the conversion from veloc- 
ity to masses. Neglecting the center-of-mass data of the 
interacting system these 14 parameters reduce to 7 pa- 
rameters containing the relative positions and velocities. 
These 7 values just fix the orbit in the case of a two-body 
interaction. Moreover, one has to specify the parameters 
that characterize both stellar systems, e.g. characteristic 
scales, orientation, or rotation. The final result is a high- 
dimensional parameter space which is in general too large 
for a standard search method. For instance, the interaction 
of a galactic disc with a point-mass galaxy is described at 
least by 7 parameters. A regular grid with a poor cover- 
age of 10 grid points per dimension demands 10'' models 
or 3400 years GRAPE simulation time (assuming 3 CPU- 
hours for a single simulation) or still about a year with a 
restricted N-body program (assuming 3 CPU-seconds per 
simulation) . 

An alternative approach to investigate large parame- 
ter spaces is a genetic algorithm (Holland 1975, Goldberg 
|1989 ). This applies an evolutionary mechanism on a pop- 
ulation (i.e. a set of parameters), from which members are 
selected for parenthood according to their fitness. Fitness 
is determined by the ability of the parameters to match 
the observations, whereas the reproduction within the ge- 
netic algorithm is done by cross-over of the parental pa- 
rameters and applying a mutation operator. Although ge- 
netic algorithms (GA) have been used in many branches of 
science, there are just a few applications in astrophysics, 
e.g. for fitting rotation curves or analysis of Doppler ve- 
locities in S Scuti stars (for a review see Charbonneau 
1995 ). Recently, Wahde ( 199^ ) demonstrated the ability 



of genetic algorithms to recover orbital parameters for ar- 
tificial data generated by N-body simulations of strongly 
interacting galaxies. He stressed that in these systems the 
positional information can be already sufficient to deter- 
mine the orbital elements. 

In this paper we combine different N-body techniques 
in order to create a model for the dynamics of interacting 
galaxies. In a first step, restricted N-body simulations are 
performed to constrain the model parameters: If a suf- 
ficiently accurate data set is available, the genetic algo- 
rithm can be used to find the orbital parameters and to 
check the uniqueness of the result. Additionally, the sen- 
sitivity of the solution can be checked by an extended 
parameter study. Alternatively, in the case of insufficient 
data, one can use the restricted N-body models to obtain 
a first guess of the parameters, creating an artificial inten- 
sity map and then apply the genetic algorithm to check 
its uniqueness. In the second step, self-consistent mod- 
els are necessary in order to tune the parameters and to 
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check the applicability of the restricted N-body method. 
These checks must address two questions: Is the neglec- 
tion of self-gravity acceptable? Does gas dynamics alter 
the results significantly? Technically this means that self- 
consistent simulations with and without gas have to be 
performed. 

The numerical models of this paper are performed 
on NGC 4449 which is a prototype actively star-forming 
galaxy similar in size and mass to the Large Magellanic 
Cloud. NGC 4449 is surrounded by an extended HI halo 
which shows some distortion (van Woerden et al. |1975; 



Bajaja et al. [1994D . Hunter et al. ( |1998D published high 
resolution HI images showing a large scale lopsided struc- 
ture. Though there seems to be a close dwarf companion, 
DDO 125, it is not clear whether the observed stream- 
ers are tidal tails caused by an interaction with DDO 125 
(Hunter et al. 199^ ). For example, tidal distortions are 
missing in the optical image of both galaxie s and no clear 
bridge has been detected. Theis & Kohle ( 1998 ) demon- 
strated qualitatively by N-body simulations that the main 
HI features could stem from a face-on encounter provided 
the mass of DDO 125 exceeds 10% of the mass of NGC 
4449. An alternative scenario assumes that the HI struc- 
ture is an ongoing infall of gas onto NGC 4449, possibly 
caused by the encounter with DDO 125 followed by a com- 
pression and cooling of the gas (Silk et al. 1987 ). However, 
it is puzzling that the observed large-scale, lopsided and 
regular structure in NGC 4449 can be maintained over 
timescales such as the orbital period which is longer than 
1 Gyr for the outer streamer. The high internal velocity 
dispersion of 10 kms~^ would destroy the streamers on 
a timescale of 2 • 10^ yr (Hunter et al. 1998). In this pa- 
per we investigate this interaction scenario, addressing the 
question of whether the observed structures can actually 
be formed by tidal interaction and which constraints on 
the mass distribution or the orbits of both galaxies can be 
derived. 

Sect. H summarizes the observational data of NGC 
4449 and DDO 125 relevant for the models of this pa- 
per. The different modeling techniques and the numeri- 
cal results are described in Sect. || (the restricted N-body 
method is introduced (Se ct. |3.l| ) and a set of first results 
is shown in Sect. 3.2 and |3.3[). The genetic algorithm ap- 



proach is briefly introduced in Sect. and explained 
in more detail in the Appendix. Its results are shown in 



Sect. 3.5. Finally, the restricted N-body simulations are 



compared with self-consistent models with (Sect. 3.7) and 
without gas (Sect. |3^ ). The results of these different ap- 
proaches are summarized and discussed in Sect. ^. 

2. The case of NGC 4449 and DDO 125 

2.1. NGC 4449 

When optical images of NGC 4449 and the Large 
Magellanic Cloud are compared, many similarities like the 
presence of a bar, the size or the total luminosity are 
found. However, unlike the LMC NGC 4449 is a candidate 



Fig. 1. Contour map of the integrated Hl-intensity of 
NGC 4449- NGC 4449 is at the center of the image and 
DDO 125 is located at the almost circular contours in the 
South. Details are given in Hunter et al. (199i ). This im- 
age was kindly provided by Deidre Hunter. 



for studies of the properties of a Magellanic-type irregular 
unaffected by the tidal field of a large galaxy. Its distance 
of 2.9 - 5.6 Mpc|^ allows for detailed investigations, per- 
formed in the visual (e.g. Crillon & Monnet 1969; Hunter 
& Gallagher |1997| ), in the radio band (e.g. van Woerden 
et al. 1975; Klein et al. 1996), in the HI emission line (e.g. 
Bajaja et al. 1994; Hunter & Gallagher [1997 ; Hunter et 



al. 1998), in the CO emission line (e.g. Sasaki ct al. 

Kohle et al. 



Hunter & Thronson |1996| Ko hle et al. |199§), in the in- 
frared (e.g. Hunter et al. |1986| ), in the FUV (e.g. HiU et al. 



1990 



1994 ) a nd in X-rays (e.g. della Ceca et al. 1997 ; Romans 
& Chu 1997). Inside the Holmberg diameter of 11 kpc, 
NGC 4449 shows ongoing s trong star formation along a 
bar- like structure (Hill et al. 1994). This activity seems to 
induce a variety of morphological features like filaments, 
arcs, or lo ops w hich have a size of several kpc (Hunter & 
Gallagher [l997| ). 

Early radio observations of van Woerden et al. ( 1975| ) 
exhibited a large HI structure exceeding the Holmberg 
diameter of 11 kp c by a factor of 7.5. Recent observations 
by Bajaja et al. ([1994 ) with the Effelsberg telescope and 
by Hunter et al. ( 1998| ) with the VLA revealed a complex 
structure of several components (Fig. mj: An elongated 



Most of the distance estimates for NGC 4449 differ be- 
cause of the assumed value Ho for the Hubble const ant. 
Independently of Ho, Karachentsev & Drozdovsky (1998) de- 



rived from the photometry of the brightest blue stars a distance 
of 2.9 Mpc. However, the error of this in easurement is unclear. 
In accordance with Hunter et al. ( |l998| ) a distance of 3.9 Mpc 
will be assumed in this paper. Thus, 1' corresponds to 1.1 kpc. 



4 



Ch. Theis & S. Kohle: Multi-Method-Modeling of Interacting Galaxies. I. A Unique Scenario for NGC 4449? 



ellipse of HI gas centered on NGC 4449 has a total mass 
of ^ 1.1 • 10^ Mq. This gas rotates rigidly inside 11 kpc, 
reaching a level of 97 kms~^ in the deprojected velocity. 
Outside 11 kpc, the rotation curve is flat. Bajaja et al. 



(|1994D derived a dynamical mass of 7.0 • 10^° Mq inside 



32 kpc and a HI gas mass fraction of 3.3%. Within 11 
kpc the dynamical mass is 2.3 • 10^° Mq and the mass 
inside the optical radius of about 5.6 kpc is 3 • 10^ Mq. 
Applyi ng sin gle dish data for a 21' beamwidth, Hunter 



et al. (|1998D derived about IO^^Mq. The HI gas inside 



the optical part of NGC 4449 counterrotates with respect 
t o the outer regions. For the central ellipse, Hunter et al. 
( 1998 ) determined a position angle of 230° ±17° and an 
inclination of 60° ±5°. 

South of the galactic center, a streamlike structure of 
25 kpc emanates which abruptly splits into two parts: A 
small spur (5 kpc) which points towards DDO 125, a close 
irregular galaxy, and a long extended stream. The latter 
first points 24 kpc to the north (in the following: vertical 
streamer), then 49 kpc north-east and, finally, 27 kpc to 
the south-east (together called the northeastern streamer), 
covering 180° around NGC 4449. The long straight sec- 
tions are notable, as are the rather abrupt changes of di- 
rection. The total mass of the streamers measured by the 



VLA observations is about 1 • 10^ Mg. Hunter et al. (1998) 



point out that compared to the Effelsberg data there is an 
additional diffuse HI mass not detected in the VLA ob- 
servations which amounts to two-thirds of the HI in the 
extended structures. 



2.2. DDO 125 

DDO 125 is located in the south-east of NGC 4449 at 
a projected distance of 41 kpc. The line-of-sight velocity 
difference between DDO 125 a nd NGC 4449 is only about 
10 kms~^. TuUy et al. (1978) found rigid rotation inside 
the optical radius of DDO 125 and an inclination of the 
adopted disc of 50°. From this they derived a total mass of 
5 • 10® Mq (corrected for a distance of 3.9 Mpc) which is in 
agreement with single dish observations (21' beamwidth) 
of Fisher & TuUy ( |1981| ). However, Ebneter et al. ( [1987D 
reported VLA-observations which show a linear increase 
of the rotation curve out to the largest radii where HI 
has been detected, i.e. out to twice the optical radius. 
This already allows for an estimate of a lower mass limit 
for DDO 125 which is a factor 8 larger than the value of 



Tully et al. (1978), i.e. 4- lO^Mg. Since neither flattening 



of the rotation curve nor any decline has been observed, 
this value is only a lower limit. If one compares that mass 
wit h the total mass of NGC 4449 derived by Bajaja et 
al. ( |l994| ) one gets a mass ratio of about q « 5.7% which 
is, however, just a lower limit. Comparison of the masses 
inside the range of rigid rotation gives a higher mass ratio 
ofq^ 18%. 



3. Numerical Modeling 

For the simulations described in this paper three different 
N-body approaches have been applied: The method of re- 



stricted N-body simulations (Pfleiderer & Siedcntopf 1961 
Toomre & Toomre 1972) is described in Sect. 3.1. The self- 



gravitating models are split into calculations without gas 



(Sect. 3.6) and simulations which include gas by means of 
sticky particles (Sect. ^ ). 

The advantage of the first method is twofold: the re- 
duction of the full N-body problem to N single body prob- 
lems and the use of particles as test particles which gives 
high spatial resolution at the places of interest for almost 
no computational cost. As a direct consequence, restricted 
N-body calculations can be performed with fewer parti- 
cles than self-consistent simulations which strongly suffer 
from two-body relaxation in case of small particle num- 
bers. Furthermore, the CPU-time scales linearily with par- 
ticle number which is superior to almost all N-body codes. 
E.g. a restricted N-body simulation sufficiently reproduces 
tidal features with only 1000 particles, whereas a direct 
self-consistent calculation needs an (optimistic) minimum 
of 10 000 particles to reproduce a stable galactic disc on 
the interaction timescale. Comparing the CPU-times of 
such a self-consistent simulation performed on a GRAPE3 
board with the values of a restricted N-body simulation we 
get a time saving of at least a factor of 1000.0 Hence, the 
restricted N-body calculation is at the moment the only 
method which allows about 10^ simulations in a couple 
of CPU-hours. Such a fast computation is a prerequisite 
for all systematic searches in parameter space such as the 



genetic algorithm described in Sect. 3.4 



However, the simplified treatment of the galactic po- 
tential in restricted N-body simulations prohibits any 
feedback of the tidal interaction on the orbits of both 
galaxies, e.g. no transfer of orbital angular momentum 
into galactic spin is possible (i.e. no merging). To overcome 
this problem, the results of restricted N-body models and 
genetic algorithms must be compared with detailed self- 
consistent simulations. In a first step the stellar dynami- 
cal applicability of the restricted models is checked by a 
comparison with self-gravitating systems calculated with 
a GR APE3af special purpose computer (Sugimoto et al. 
19901; see Sect. IsT 



So far the dynamics of the system have been inves- 
tigated in terms of stellar dynamics, though in many 
cases (such as the case of NGC 4449) HI gas is observed. 
In general, the effects of neglecting gas dynamics are 
not clear. E.g. a purely stellar dynamical ansatz was 
successfully applied for the s ystem M 81 and NGC 3077 
(Thomasson & Donner 1993|) , whereas in other systems 
the dynamics of the gas can deviate s trong ly from the 
behaviour of the stars (e.g. Noguchi 1988| ; Barnes & 



^ The exact acceleration depends strongly on the details of 
the comparison, e.g. the requirements on disc stability and 
therefore the number of particles for the self-consistent model, 
the detailed orbits, the method of the self-consistent calcula- 
tion etc. An acceleration of a factor 1000 is only a lower limit. 
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Hernquist 1996). Therefore, one has to compare stellar 
and gas dynamical results which is done in this paper by 
a sticky particle method (Sect. 3.7). 



The units are chosen as 1 kpc, 7 • 10 Mq (the dynam- 
ical mass of NGC 4449) and G=l. This results in a time 
unit of 1.78 Myr and a velocity of 549 kms~^. 



3.1. Restricted N-body models 

Restricted N-body simulations aim to reduce the N-body 
problem to N 1-body problems by assuming that the grav- 
itational potential is given by a simple relation. E.g. two- 
or a few-body problems have known (semi-)analytical so- 
lutions or can be solved by fast standard methods. Here 
we assume that the gravitational forces on each particle 
are given by the superposition of the forces exerted by two 
point-like objects (the galaxies) moving on Keplerian or- 
bits. Thus, for the orbits one has to specify 14 parameters 
which reduce to 7 in the center-of-mass frame (relative 
coordinates/ velocities and the mass ratio): The orbital 
plane is fixed by the inclination angle and the argument 
of the ascending node. The orbit itself is characterized by 
its eccentricity and the location of the pericenter, i.e. the 
pericentric distance and the argument of the pericenter. 
Finally, mass and initial location of the reduced particle 
(or the time of pericentric passage) have to be specified to 
fix the phase. For a comparison with the observations, the 
phase of the observation also must be supplied. Three of 
the seven orbital parameters can be fixed by the observed 
location of both galaxies on the plane of sky {x-y plane) 
and the line-of-sight velocity. The masses (and thus the 
mass ratio) can be estimated by the velocity profiles of 
both galaxies, provided the distance to them is known. 

In this paper, the time of pericentric passage is set 
to zero. As input, the orbital eccentricity and inclination 
and the minimum distance are given. Furthermore, the di- 
rectly observable parameters, i.e. the masses of the galax- 
ies, the projection of the relative positions of both galax- 
ies onto the plane of sky, and their line-of-sight velocity 
difference, are supplied. They are used to determine the 
argument of the pericenter, the phase (or time since peri- 
centric passage) of the observed system and the argument 
of the ascending node of the orbital plane. Additionally, 
the initial distance or phase has to be specified: In the case 
of parabolic and hyperbolic encounters, the distance was 
set to 100 kpc, which guarantees a sufficiently low tidal 
influence at the start of the simulation. In the case of an 
elliptical encounter, the situation is more difficult, because 
the apocentric distance might be small, especially if the 
eccentricity is low. Thus, the system might be tidally af- 
fected all the time or by a series of repeated encounters 
which makes the applicability of the restricted N-body 
method doubtful. For these elliptical orbits, the initial az- 
imuthal angle of the orbit was set to 120° (0° corresponds 
to pericenter). 



The test particles are arranged in a flat disc moving on 
circular orbits. The disc itself is characterized by its ori- 
entation, i.e. its inclination and position angle, the scale- 
length and the outer edge. Additionally, an inner edge can 
be specified in order to avoid a waste of computational 
time by integrating tidally unaffected orbits well inside 
the tidal radius of the galaxy. In most of the restricted 
N-body simulations the test particles are distributed in 
30 to 60 equidistant rings of 125 particles, i.e. a total of 
N — 3750 to 7500 particles. The typical spacing of the 
rings is 0.7 to 1.5 kpc depending on the outer edge of 
the disc. The time integration of the equations of motion 
is performed by a Cash-Karp Runge-Kutta integrator of 
fifth order with adaptive timestep control (Press et al. 
1992). In order to prevent a vanishingly small timestep by 



an accidental infall of a test particle to the galactic center, 
the gravitational force is softened on a length scale of 1 
kpc. In order to speed-up the calculation the galactic po- 
sitions are tabulated by cubic splines and the integrator 
uses these look-up tables. On a SPARCIO with a 150 MHz 
CPU a typical 7500-particle encounter takes 23 CPU sec- 
onds, which can be reduced by a factor of 2-3 by applying 
an inner edge of 10 kpc. 

3.2. A first guess 

In a series of simulations the orbital parameters (minimum 
distance d„iin, the eccentricity e, the mass ratio q and the 
orientation of the orbital plane) as well as the disc pa- 
rameters (size i?max of the disc, its position angle a and 
inclination i) are varied (see Table [sT^ ) . 

The evolution of the primary galaxy (NGC 4449) in the 
reference model A is displayed in Fig. |^. At the beginning, 
both galaxies have a distance of 100 kpc corresponding to 
995 Myr before closest approach. The particles inside 10 
kpc are omitted for an easier identification of the evolving 
structure. This does not alter the results, since the related 
particles are only weakly affected by the interaction, as 
the persistence of the elliptical shape of the central hole 
demonstrates. At the moment of closest approach {t = 0), 
only the outer region of NGC 4449 's disc reacts on the 
intruder. 71 Myr later (i.e. system time 40), a dense outer 
rim, starting at the projected position of DDO 125 and 
pointing to the north, begins to form. Secondly, a weak 
bridge-like feature seems to connect both galaxies. Both 
features are more pronounced after 178 Myr. At that time 
an additional tidal feature begins to evolve north of the 
center. When DDO 125 has a projected distance of 41 
kpc (362 Myr after closest approach), three features are 
discernible: the remnant of the former bridge starts south 
of the center of NGC 4449 pointing to the south-west. 
The vertical feature starts at the end of the first structure 
and points to the north. The third structure is a long tidal 
arm which starts west of NGC 4449 pointing to north-east 
before turning back to the south-east. 

A comparison of the intensity map of model A (lower 
right diagram in Fig. ||) with the observations (Fig. |lj) 
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Fig. 2. Temporal evolution of the particle positions around 
NGC 4449 projected onto the plane-of-sky for different 
times. The perturber DDO 125 is shown by the star. 
Closest approach is at t — 0. The final diagram corre- 
sponds to a projected distance of 4I kpc. The mass ratio 
of both galaxies is set to q ^ 0.2. The original HI disc 
has an inclination of i — 60°, a position angle a — 230° 
and a radial extent of i?max = 40 kpc. The particles in- 
side a radius of 10 kpc are omitted. The contour lines in 
the lower right diagram show the normalized surface den- 
sity distribution corresponding to the particle distribution 
in the lower left diagram. The intensities are calculated 
on a 50x 50 grid and .smoothed over 3 grid cells in each 
dimension. One time unit corresponds to 1.78 Myr. 



demonstrates that the main observational features can be 
reproduced qualitatively by the reference model A, i.e. a 
parabolic orbit, a mass ratio q = 0.2 and a minimum dis- 
tance dmin = 25 kpc. The parameters for the disc orienta- 
tion are chosen in agreement with Hunter et al.'s sugges- 
tion oi i — 60°, a = 230°, whereas the disc radius was set 
to i?max = 40 kpc. Although this numerical model is not 
a detailed fit of the data, the characteristic sizes and lo- 
cations of the streamers as well as the unaffected disc can 
be found. The material located between the streamer and 
the disc, which is not seen in Hunter's data, corresponds 
to the extended HI mass detecte d in th e single dish ob- 
servations by van Woerden et al. ( 1975|) . In the following 
sections we will discuss the influence of the individual pa- 



Fig. 3. Dependence on the eccentricity. Projection of the 
particle positions onto the plane-of-sky at the moment of 
a projected distance of 4I kpc. The eccentricities of the 
models are 0.5 (upper left, model Bl), 0.8 (upper right, 
B2), 1.0 (A, reference model), 1.5 (33) and 5.0 (B4). 



rameters. The global uniqueness of the reference model is 



investigated in Sect. 3.5 



3.3. Parameter study 

In this section the influence of the basic parameters, 
i.e. the orientation of the orbital plane, the orbital 
eccentricity and the minimum distance, the mass ratio 
of both galaxies, and the orientation of the HI disc and 
its extent are discussed (see Table 3.3). In the following 
simulations, all parameters except the one of interest are 
kept constant with respect to model A. Thus, the vicinity 
of the reference model in parameter space is studied. 

Eccentricity. The eccentricity mainly influences the in- 
teraction timescale which is defined here as the time be- 
tween maximum approach and a projected distance of 41 
kpc, i.e. the moment of observation. An eccentricity of 5.0 
(corresponding to a hyperbolic encounter) decreases the 
interaction time by a factor 2.6 resulting in less structural 
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Model parameters of the restricted N-body simulations 



model 


e 


<? 


drain 


i 


id 




^max 


comment 


A 


1.0 


0.2 


25 


40 


60 


230 


40 


reference model 


Bl 


0.5 


0.2 


25 


40 


60 


230 


40 


eccentricity 


B2 


0.8 


0.2 


25 


40 


60 


230 


40 




B3 


1.5 


0.2 


25 


40 


60 


230 


40 




B4 


5.0 


0.2 


25 


40 


60 


230 


40 




CI 


1.0 


0.01 


25 


40 


60 


230 


40 


mass ratio 


C2 


1.0 


0.05 


25 


40 


60 


230 


40 




C3 


1.0 


0.1 


25 


40 


60 


230 


40 




C4 


1.0 


0.3 


25 


40 


60 


230 


40 




C5 


1.0 


0.4 


25 


40 


60 


230 


40 




Dl 


1.0 


0.2 


15 


40 


60 


230 


40 


minimum distance 


D2 


1.0 


0.2 


20 


40 


60 


230 


40 




D3 


1.0 


0.2 


30 


40 


60 


230 


40 




D4 


1.0 


0.2 


35 


40 


60 


230 


40 




D5 


1.0 


0.2 


40 


40 


60 


230 


40 




El 


1.0 


0.2 


25 


10 


60 


230 


40 


orbital inclination 


E2 


1.0 


0.2 


25 


20 


60 


230 


40 




E3 


1.0 


0.2 


25 


30 


60 


230 


40 




E4 


1.0 


0.2 


25 


50 


60 


230 


40 




E5 


1.0 


0.2 


25 


70 


60 


230 


40 




Fl 


1.0 


0.2 


25 


40 


40 


230 


40 


disc inclination 


F2 


1.0 


0.2 


25 


40 


50 


230 


40 




F3 


1.0 


0.2 


25 


40 


70 


230 


40 




Gl 


1.0 


0.2 


25 


40 


60 


210 


40 


pos. angle of disc 


G2 


1.0 


0.2 


25 


40 


60 


250 


40 




HI 


1.0 


0.2 


25 


40 


60 


230 


25 


outer edge of disc 


H2 


1.0 


0.2 


25 


40 


60 


230 


30 




H3 


1.0 


0.2 


25 


40 


60 


230 


35 




H4 


1.0 


0.2 


25 


40 


60 


230 


50 





Table 1. Parameters of the different restricted N-body models: model name (col. 1), orbital eccentricity e (col. 2), 
mass ratio q of both galaxies (col. 3), minimum separation dmin of galactic centers (in kpc, col. 4), inclination i of 
orbital plane (in degrees, col. 5), orientation of NGC 4449 disc (in degrees): inclination id (col. 6) and position angle 
ad (col. 1), maximum radius i?max of gaseous disc (in kpc, col. 8). 



changes of the primary galaxy (Fig. |^): Only the rem- 
nant of the bridge feature can be discerned, whereas the 
other two streamers found in the reference model are miss- 
ing completely. With decreasing eccentricity, the north- 
eastern streamer becomes pronounced before the vertical 
structure arises. The latter is weakly expressed for e = 1.5. 
Reducing e further to bound orbits makes the features 
even more prominent. In the case of e = 0.5 the north- 
eastern and the vertical streamer form a single structure 
clearly detached from the center. However, due to the 
strength of the interaction, a second broad tidal arm is 
ejected from NGC 4449, leaving a low-density region south 
to the northern tip of the north-eastern streamer. 

In order to compare the e = 0.5 model with observa- 
tions, we reran model Bl, but with = 80 000, in order to 
derive intensity and velocity maps. This model reproduces 
several observed features very well (Fig. U): First, the in- 
tensity map shows the three observed abrupt changes: at 
the lower tip of the vertical streamer, at the transition 
from the vertical streamer to the north-eastern streamer 



and at the northern tip of the north-eastern streamer. 
Second, the spatial extent of the streamers is in agreement 
with the observations. Third, we find a bridge between the 
north-eastern streamer and the central tidally-unaffected 
part of NGC 4449. Fourth, during the encounter, mate- 
rial is ejected south to NGC 4449. The latter two features 
can be found in deeper images of NGC 4449 (cf. Fig. 1 of 
Hunter et al. 19981 ). However, there are also some remark- 
able differences: First, both model and observation show 
a spur at the lower tip of the vertical streamer. However, 
their directions differ by about 90°. Second, the position 
of the northern tip of the north-eastern streamer (or the 
angle between vertical and the north-eastern streamer) 
differ. 

Looking at the velocity map, we find Keplerian 
rotation in the central region of the simulations, whereas 
the observations show rigid rotation. The velocities 
in the streamers reach about 80-100 kms""'^, which is 
in agreement with the observed values. Additionally, 
the velocity field in the simulations shows a small. 
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Fig. 4. Projected particle positions (upper left), velocity 
map (upper right) and intensity map (lower left) for model 
Bl (e — 0.5) at the moment of a projected distance of 41 
kpc. The velocity contours differ by 20 kms^^, whereas the 
intensity map ranges from 1% of the maximum intensity 
Imax to /max- The star marks the position of DDO 125. 
The observed intensity map ( cf. Fig. 1) is displayed in the 
lower right panel for comparison. 



but significant V-shape, south of NGC 44 49 wh ich is 
also observed (cf. Fig. 3 of Hunter et al. |l998D. The 



main discrepancies between the observed and numerical 
velocity maps comes from the point mass approximation 
of the galactic potential used for the numerical model. 
Although this approximation must fail for the central 
region, it should be sufficient for the outer regions which 
are most affected by tides and which contain almost the 
entire dynamical mass of NGC 4449. 

Mass ratio. The variation of the mass ratio gives a lower 
limit of about 10% in order to create a sufficient tidal 
response (Fig. ||). If the mass ratio exceeds 0.3, the ver- 
tical feature becomes less vertical and the north-eastern 
streamer is more diffuse, both due to the enhanced tidal 
pull of the intruder. 

Minimum distance. Although the interaction time 
is not strongly influenced by the minimum distance 
drain, it determines the maximum strength of the tidal 
field (Fig. ||): For small separations, the structures 
become more diffuse and the primary is much more 
affected. The vertical feature is more extended, creating 
an additional tidal arm (similar to the model with an 
orbital eccentricity e — 0.5). If the distance of closest 
approach exceeds 30 kpc, almost no streamers are formed. 



Fig. 5. Dependence on the mass ratio. Projection of the 
particle positions onto the plane-of-sky at the moment of a 
projected distance of 41 kpc. The mass ratios are - starting 
with the upper left diagram - 0.01 (model CI), 0.05 (C2), 
0.1 (C3), 0.2 (A, reference model), 0.3 (C4) and 0.4 (C5). 



Orbital plane. The orientation of the orbital plane char- 
acterizes the observer's location or the line-of-sight, which 
is fixed by two parameters. For the models here only the 
orbital inclination is a free parameter, since the second pa- 
rameter, the position angle of the orbital plane, is fixed by 
the observed relative line-of-sight velocity of both galaxies. 
The models show that the inclination affects mainly the 
orientation of the vertical streamer (Fig. |^): An inclina- 
tion outside the range [30°, 50°] gives too much deviation 
from the south-north direction of the vertical feature. 

HI distribution in NGC 4449. A specification of the 
mass- and especially the HI distribution is a more diffi- 
cult task, because, contrary to the well-known two-body 
problem, the number of free parameters characterizing the 
structure and internal dynamics of the interacting galaxies 
is unknown and a matter of investigation itself. For exam- 
ple, it is not a priori clear how the HI is distributed. To 
reduce the number of free parameters, we will assume here 
that the observed inner elongated ellipse is the remnant 
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Fig. 6. Dependence on the minimum distance during the 
encounter. Projection of the particle positions onto the 
plane- of-sky at the moment of a projected distance of 41 
kpc. The minimum distances are - starting with the upper 
left diagram - 15 kpc (model Dl), 20 kpc (D2), 25 kpc (A, 
reference model), 30 kpc (D3), 35 kpc (D^) and 40 kpc 
(D5). 



Fig. 7. Dependence on the inclination angle of the orbital 
plane. Projection of the particle positions onto the plane- 
of-sky at the moment of a projected distance of 41 kpc. The 
inclinations are Iff' (upper left, model El), 20^ (upper 
right, E2), 30^ (E3), ^0° (A, reference model), 50° (E4) 
and 70° (E5). 



of an original exponential HI disc, which can be specified 
by its orientation, its scale length and radial extent. The 
inclination of the disc can roughly be determined from 
the HI map, whereas its sign is not unambiguous. Placing 
the south-eastern side of the HI disc into the foreground 
means that the orbital plane of the companion and the 
plane of the HI disc are nearly orthogonal. We performed 
several simulations with geometries of this kind, but these 
show almost no strong tidal responses and they all result 
in situations completely different to that observed (Kohle 



199£). However, placing the north-western side of the HI 



disc into the foreground and considering the HI velocity 
field, we now have the companion on a prograde orbit. 

Since the particles are treated as test particles, a vari- 
ation of the scale length does not change the final location 
of a particle, but rather its mass, which is required for the 
calculation of the intensity map. A small scale length (e.g. 
5 kpc) would confine too much HI in the center which is 



unaffected by the encounter. Since the observations show 
a large fraction of HI outside 20 kpc, we have chosen a 
large scale length of 30 kpc which gives only a weak radial 
decline of the HI surface density and, by this, a sufficient 
amount of HI in the outer regions. 

Disc size. The influence of the maximum radius i?max of 
the disc is shown in Fig. |^. If i?inax is smaller than 35 kpc, 
the vertical streamer has not been formed at all, indicating 
that this feature stems from the HI gas far outside the 
center. The north-eastern arm and the bridge, however, 
are already built up for i?max = 25 kpc. With increasing 
disc size the vertical streamer becomes more pronounced 
and more extended, which excludes disc sizes larger than 
45 kpc. 

Disc orientation. The orientation of the HI disc also 
strongly affects the appearance of the primary galaxy. 
Variation of the position angle a by 20° gives HI dis- 
tributions which are incompatible with the observations 
(Fig. ^. If a is decreased, the vertical streamer is more 
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Fig. 8. Dependence on the extension of the HI distribution 
in NGC 4449- Projection of the particle positions onto the 
plane- of- sky at the moment of a projected distance of 41 
kpc. The radial extensions are 25 kpc (upper left, model 
HI), 30 kpc (upper right, H2), 35 kpc (H3), 40 kpc (A, 
reference model) and 50 kpc (H4). 

extended and diffuse. On tlie otlier hand, if a is increased, 

the north-eastern arm becomes too long and the vertical 
feature vanishes almost completely. Similarly, the results 
depend strongly on the inclination i: Inclinations below 
50° give only a weak bridge and a weak vertical streamer, 
whereas the orientation of the vertical streamer deviates 
strongly from the northern direction for inclinations 
exceeding 70°. 

The performed parameter study demonstrates that 
even close to the parameters of model A, the final par- 
ticle distribution shows significant variations which are 
not in agreement with the observations. Thus, the ref- 
erence model might be a good and locally unique candi- 
date describing the dynamics of NGC 4449 and DDO 125. 
The confidence in this encounter scenario is especially en- 
hanced by the strong dependence on the orientation of the 
HI disc of NGC 4449. Qualitatively good models are only 
found if the disc orientation agrees well with the values 



Fig. 9. Dependence on the orientation of the HI disc in 
NGC 4449- Projection of the particle positions onto the 
plane- of- sky at the moment of a projected distance of 41 
kpc. The orientation is characterized by the inclination id 
of the disc and the position angle ad'- id = 60°, ad = 230° 
(upper left, model A, reference model); 6(f , 21 CP (upper 
right, Gl); 6(P , 25ff' (G2);4(P, 23CP (Fl); 5CP , 23(P (F2) 
and 70°, 23CP (F3). 

determined independently by observations. On the other 

hand, the former parameter study only investigates a small 
fraction of parameter space and, thus, the uniqueness of 
the solution is not confirmed. Moreover, in general, it is 
very tedious to find a good model by hand, i.e. by an 
unsystematic non-automated search in parameter space. 
Thus, the probability of accepting the first solution found 
seems to be very high (especially, if CPU-intensive simu- 
lations are performed), by this neglecting other regions in 
parameter space which might also give sufficient fits. In the 
next section an automatic search in parameter space by 
means of a genetic algorithm is introduced. This allows in 
principle for an ab initio determination of acceptable pa- 
rameter sets (if data sets of sufficient quality are available) 
or at least for a uniqueness test of a favoured scenario like 
the reference model described above. 
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3.4. The genetic algorithm 

The idea of applying models of organic evolution for opti- 
mization pro blems dates bac k to th e 1960s and 1970s (e.g. 
Rechenberg 1965; Schwefel 1977). Unlike standard de- 



terministic gradient techniques for optimization (e.g. the 
downhill simplex method (Press et al. 1992[ )) Rechenberg's 
Evolutionsstrategie is probabilistic: Starting with a more 
or less random parent^ i.e. a single point in parameter 
space, a child is generated by a random mutation of the pa- 
rameter set characterizing the parent. The quality of both 
individuals with respect to the optimization problem (i.e. 
their fitness) determines the parent of the next generation. 
Repeating this process of mutation and selection improves 
the quality of the individual monotonously. This special 
implementation of evolutionary algorithms (EAs) in gen- 
eral, has certain advantages: Compared to complete grids 
in parameter space, the probabilistic but oriented nature 
of the evolutionary search strategy allows for an efhcient 
check of a high-dimensional parameter space. Compared 
to gradient methods which are very fast near the optimum, 
EAs do not need any gradient information which might be 
computationally expensive. They depend only weakly on 
the starting point and - most important - they are able to 
leave local optima. The price for these features is a large 
number of fitness evaluations or test points in parameter 
space before converging to a good solution. 



Fig. 10. Best fit model during the course of a GA fitting 
procedure. The projection of the particles in x-y-plane and 
the corresponding grid for the intensity evaluation are dis- 
played: The original data (upper left), the best fit of the GA 
after initialization ( upper middle ), after the first breeding 
(upper right), after 11 generations (lower left) and at the 
end of the fitting procedure after 100 generations (lower 
middle). The evolution of the maximum fitness is shown 
in the lower right diagram. The number of test particles 
per simulation is 900. 



Holland (197f;) extended the Evolutionsstrategie of Re- 
chenberg by abandoning asexual reproduction: He used 
a population of individuals instead of a single individual 
(and its offspring). From the population two individuals 
were selected according to their fitness to become par- 
ents. These parents represent two points {chromosomes) 
in parameter space and the corresponding coordinates are 
treated like genes on a chromosome. These chromosomes 
are subject to a cross-over and a mutation operation re- 
sulting in a new individual which is a member of the next 
generation. Such a breeding is repeated until the next gen- 
eration has been formed. Finally, the whole process of sex- 
ual reproduction is repeated iteratively until the individu- 
als representing a set of points in parameter space confine 
sufficiently one or several regions of high fitness. Despite 
many differences in the detailed implementations of ge- 
netic algorithms (GAs), the basic concept of an iterated 
application of randomized processes like recombination, 
mutation, and selection on a population of individuals re- 
mains in all GAs (Goldberg |1989| ; Back |1996D . 



3.5. Are the results unique? 

In order to check the uniqueness of the reference model 
A, several runs with different parameter sets encoded in 
genes were performed. Fig. ^ shows a typical run of the 
GA operating on a population of 100 members. The ref- 
erence model^ - displayed in the upper left diagram - as 
well as the GA models were calculated with 900 test par- 
ticles. This makes the streamers less prominent, but still 
discernible. The next plots show the GA status at different 
generations: At generation 1, the best model is created out 
of the randomly chosen initial parameter sets. However, 
already by the second generation, i.e. after the first ap- 
plication of the reproduction operator, the best model 
clearly exhibits the north-eastern streamer. At generation 
11, even the weaker vertical feature and the bridge feature 
appear. At the end of this GA run the particle distribu- 
tion of the original model A and the best fit of the GA are 
almost identical, at least compared by eye. The fitness of 
the best fit has strongly increased from 0.04 to 0.14 where 
values above 0.1 indicate a very good fit. Most of the im- 
provements were found during the first 50 generations, i.e. 
analyzing 5000 models. After that generation, the homo- 
geneity of the population led to inbreeding which strongly 
prohibits the appearance and spread of new (better) pa- 
rameters. However, this is not a problem for the models 
here, because a very good fit to all parameters was already 
found after 40 generations (Fig. |ll|). The relative devia- 
tion of the derived parameters from the original values is 
less than 15% in all cases, and for many is better than 5%. 

In a series of different GA runs, we varied the pop- 
ulation size, the number of generations or particles, the 



We were only able to perform a uniqueness check of our 
preferred scenario using the best model as reference. A better 
model might be determined using the original VLA data as 
reference. 
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Fig. 12. Test of the capability of the GA (lower row) to 
recover original data when these data (upper row) are of 
low quality or incomplete. Shown are runs with low reso- 
lution, i.e. a discretization on a 4^4 QT^d (left column), a 
had pixel (here an omitted central pixel, middle column), 
and a missing structure (the vertical streamer, right col- 
umn). 



Fig. 11. Development of six parameters of the best fit 
model during a GA run. The parameters are the mass of 
the secondary galaxy (upper left), the minimum distance 
(upper right), the inclination (middle left) and position 
angle (middle right) of the orbital plane, and the inclina- 
tion (lower left) and position angle (lower right) of the 
disc of the primary galaxy. The filled squares show the 
parameters of the fitted artificial reference model. 



parameter combination encoded on ciiromosomes, tlie fit- 
ness function, and the random initial population. In al- 
most all of them the best fit gave an acceptable fit of 
the reference model which demonstrates the capability of 
GAs to recover the parameters, even for weak encounters. 
The models which failed suffered from inbreeding or from 
a fitness function which did not discriminate sufficiently 
between high and low quality fits. Moreover, a different re- 
gion in parameter space producing a good fit (/ > 0.1) was 
never found. Thus, the model A - or more exact the cor- 
responding intensity map - seems to result from a unique 
parameter set. However, this is only a motivated guess, 
not a mathematical proof! 

Another interesting question is how strongly do the 
GA results depend on the completeness or accuracy of the 
data. We studied this in two ways: First, we used a coarser 
grid (4x4) in order to mimic a lower resolution of obser- 
vational data (left column in Fig. [ij). The GA is then 
only qualitatively able to reproduce the reference model. 
Obviously, too much information is lost when the resolu- 



tion becomes too small. On the other hand, an increase of 
the number of cells to 10x10 did not improve the fit found 
on a 7x7 grid, but introduced more noise into the intensity 
maps. In the second test we discarded one or several grid 
cells from the fitness calculation (middle and right column 
in Fig. |l^. Omitting one cell does not affect the final fit, 
independent of the location of the cell (here, the center). 
However, if a whole feature, such as the vertical streamer, 
is neglected, the GA is unable to recover the original data. 

3.6. Self-consistent models 'sin gas' 

The self-consistent models in this paper were performed by 
a direct N-body integration using the GRAPESaf board in 
Kiel for the calculation of the gravitational forces and the 
potentials. The time integration was done by a leap-frog 
scheme with a fixed timestep of 1.78 Myr. The GRAPE 
board uses intrinsically a Plummer softening, the soften- 
ing length was chosen to be 0.2 kpc. These choices give a 
total energy conservation of better than 1% for the whole 
simulation. 

Different to the restricted N-body models, the set-up 
of the initial particle configuration for self-gravitating sys- 
tems is far from being trivial. In principle, any stable sta- 
tionary stellar system fulfills the coUisionless Boltzmann 
equation. Thus, a solution of the latter could be used as 
a trial galactic system. From observations of the Milky 
Way it is known, however, that even if the Galaxy is ax- 
isymmetric, three integrals of motion are required t o con - 
struct a stationary model (e.g. Binney & Tremaine 1987 ). 
Unfortunately, only two are known analytically (the en- 
ergy and the z-component of the angular momentum). 
Additionally, even if a particle configuration is generated 
from a distribution function depending on the integrals of 
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Fig. 13. Comparison of restricted N-body models (left) 
with self- consistent simulations (right) for a parabolic en- 
counter (upper row) and an eccentricity of e = 0.5 (mid- 
dle row). The projected distance to the perturber (here 
about 40 kpc) defines the time of the snapshot. The lower 
row displays the normalized intensity profiles for the self- 
consistent simulations (left: e = 1, right: e — 0.5). Note 
that the surface plots are seen from north-west in order to 
emphasize the coherent streamer structure. 



motion, it is unclear whether it is stable against pertur- 
bations which might lead to e.g. the Toomre instability or 
the bar instability. Thus, on the timescale of interest one 
has to check the stability by numerical simulations. 

In this paper we use the models of Kuijken & Dubinski 
( |l995| ) consisting of three components which are deter- 
mined self-consistently from the Boltzmann equation: the 
bulge follows a King model and the halo a lowered Evans 
model. For the disc the third integral is approximated by 
the z-energy, i.e. the energy in vertical oscillations. The 
advantage of Kuijken & Dubinski's method is the possibil- 
ity to specify parameters hke scale length or scale height 
of the disc directly, whereas many other methods (e.g. 
Barnes 1988) not starting in exact virial equilibrium read- 
just the mass distribution on a dynamical timescale, which 
means less control on the structural parameters. 

Our initial model of NGC 4449 is derived from the ob- 



servational data of Bajaja et al. (1994) and Hunter et al. 
( |1998D . The total mass is set according to the dynamical 
mass of 7 • 10^" Mq within 32 kpc. 90% of this mass is 
attributed to the dark halo, 3% to the disc and 7% to the 
bulge component. The radial extent of the halo was set 
to 32 kpc. The disc modelled here should not be confused 
with the optical disc of the galaxy. Since we are interested 



Fig. 14. Self- consistent N-body simulation 'sin gas': 
Projection of the particle positions onto the plane-of-sky at 
different times: initially (upper row), at closest approach 
(middle row) and at the projected distance of J^l kpc (lower 
row). The left column shows the disc particles, whereas the 
right column displays the halo particles. 



in the tidal response of the outer regions, especially the 
formation of the streamers, we did not model the inner 
(optical) region of NGC 4449, i.e. the stellar disc. The 
missing signatures of any tidal features in the optical im- 
age demonstrate that the stellar body of NGC 4449 is 
not affected by the interaction. For the same reason we 
did not consider the counter-rotation of the inner region. 
However, its gravity is taken into account by means of the 
bulge and halo contributions to the overall rotation curve. 
The radial extension of the bulge is 3 kpc. The scale length 
of the extended (HI) disc was varied between 10 and 30 
kpc, the outer edge was set to 40 kpc. 

The long-term stability was checked by simulations of 
an isolated disc with N = 18000 particles (8000 for the 
disc, 4000 for the bulge and 6000 for the halo). Within 
1.5 Gyr the Lagrange radii at 10% and 90% mass vary 
only by 2-4%. They show no systematic trend except for a 
small expansion during the first 100 Myr which is proba- 
bly caused by gravitational softening. Due to two-body 
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relaxation the scale height of the disc increases by al- 
most a factor of 2 within 1.5 Gyr. However, this heating is 
not expected to alter the result of the dominant interac- 
tion significantly. The Fourier amplitude C2, defined by 
Crn = (^i,,E(i?»rdre-'"^d</>)/(Xii,^E(i?,</.)rdrd</>)) 
(with the surface density E), is always less than 5%, in 
agreement with no significant bar or spiral pattern built 
up during the simulation. 

Fig. ^ shows a comparison of restricted N-body sim- 
ulations with self-consistent models. Except for the more 
diffuse structure in the latter, the large-scale features are 
very similar in both simulations. The streamer structure 
can be found in both self-consistent simulations, as the 
surface plots demonstrate. The fact that the streamers 
are better resolved in the simulation with an eccentricity 
e = 0.5, is partly due to the increased number of par- 
ticles (80 000) used for this simulation. The simulations 
for the parabolic orbits were performed with 8 000 (disk) 
particles for both, the restricted and the self-consistent 
simulations. The latter clearly demonstrates the ability of 
the restricted N-body models to resolve structures, even 
when a small number of particles is used. It is remark- 
able that the extended dark halo which has been applied 
for the self-consistent models does not strongly affect the 
results as the comparison of the e = 0.5-models shows: 
The vertical and the north-eastern streamer, the mass 
ejected to north of the north-eastern streamer, the small 
bridge between the central region and the eastern tip of 
the streamer, the extent of the HI distribution south and 
south-east to the center show up in both simulations at 
almost the same locations. The only significant difference 
is that the north-eastern streamer is less straight in the 
self-consistent model than in the restricted N-body simula- 
tion. This very good agreement confirms the applicability 
of the restricted N-body method here. Another interesting 
aspect concerns the response of the halo to the intruder. 
Unlike to the disc, the halo remains almost unchanged by 
the encounter (Fig. |lj). This behaviour allows the appli- 
cation of a rigid galactic potential which has been used in 
the restricted N-body simulations. 



3.7. Self-consistent models 'con gas' 

For the gasdynamical models we use a sticky particle 
scheme which treats each particle of the HI disc as a 
gaseous cloud. The radius of these clouds is assumed to 
follow i?ci = 20 • v'Mci/IO'^Mq pc which is de rived from 



the mass-radius relation of Rivolo & Solomon (1987) cor- 
rected for glancing collisions (Md is the mass of a gas 
cloud). When the distance of two clouds falls below the 
sum of their radii, merging of the clouds is allowed, if 
their relative orbital angular momentum is smaller than 
the maximum angular momentum a merged cloud can ac- 
quire before break-up. If this angular momentum crite- 
rion prevents a merger, the clouds keep their kinematical 
data, i.e. no energy is dissipated. Star formation is only 
qualitatively included by adopting an upper mass limit of 



Fig. 15. Projection of the particle positions onto the 
plane-of-sky at different times: at closest approach (up- 
per row), at the projected distance of 4I kpc (middle row) 
and at the end of the simulation (lower row). Both models 
include self- gravity. The left column shows a purely stellar 
simulation, whereas the right column includes gas dynam- 
ics by means of sticky particles. 



7 • 10^ M0 for clouds. If they exceed this limit no further 
merging is allowed, i.e. they are treated as stars. However, 
the star formation criterion is not important for the mod- 
els of this paper. (For details of the cloud model or the 
merging mechanism see Theis & Hensler 1993). This dis- 
sipative scheme has been successfully applied to collaps- 
ing systems reproducing many observed properties of spi- 
ral and elliptical galaxies. For the long-term evolution of 
barred disc galaxies, a comparison of this scheme with 
that of Palous et al. ( 1993| ) for their standard coefficients 
of restitution of Pr — Pt — 0.2 gives basically the same re- 
sults (Jungwiert, private communication). Thus, the cho- 
sen model should be applicable for a reliable determina- 
tion of the difference between gas dynamical and purely 
stellar dynamical models. 



A comparison of gaseous models with dissipationless 
models shows no significant difference (Fig. ^). Due to 
the low gas densities only about 20-30% of the clouds un- 
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dergo inelastic collisions which mainly occur in the denser 
central regions. Thus, the outer HI halo is only very weakly 
influenced by (sticky particle) gas dynamics and the appli- 
cability of purely stellar simulations for the HI is possible. 

4. Discussion 

4.1. Interaction between NGC 4449 and DDO 125? 

The previous simulations demonstrate that models based 
on an encounter scenario are able to reproduce the mor- 
phology of the three streamers. Furthermore, the inter- 
action puts gas to the south-west of NGC 4449, where a 
clumpy HI distribution is observed (Hunter et al. 1998| ). 
Though its detailed patchy structure has not been recov- 
ered by the models, an initial gas distribution less regular 
than the smooth disc-like one used for the numerical mod- 
els might be responsible for this deviation. 

The parameter study showed that small variations in 
the parameters give HI distributions which are not in 
agreement with the observations. Especially, the strong 
dependence of the final structure on the initial orienta- 
tion of the HI disc confines its inclination and position 
angle within 10° to the values derived independently from 
the observations. This remarkable coincidence would oc- 
cur necessarily in the encounter scenario, whereas in the 
infall scenario it happens just by chance. A second critical 
feature is the abrupt direction change of the streamers. In 
the infall scenario they mark the path of the perturbed 
gas which collapses to the galactic center. However, it is 
hard to imagine how a realistic static potential could cre- 
ate orbits with such strong direction changes as e.g. the 
one from the vertical to the bridge feature. Moreover, due 
to the large distance of these features from the galactic 
center the potential should be more or less spheroidal or 
already spherical, which makes it even more difficult to 
produce such orbits. 

The most critical parameter for the interaction sce- 
nario is the mass ratio q of DDO 125 and NGC 4449. The 
simulations demonstrate that a mass ratio smaller than 
10% cannot reproduce the streamer structure. In order to 
determine q, both galactic masses have to be known. In 
Bajaja et al.'s ( |1994[ ) estimate for the mass of NGC 4449 
they assumed that the HI is in rotational equilibrium even 
at a distance of 32 kpc. However, in the case of an inter- 
action with DDO 125 the velocities at large distances are 
strongly affected by the interaction and they are definitely 
not in rotational equilibrium (as the ejection of the tidal 
arms in a continuation of the simulations demonstrates). 
Hence, the derived mass of 7 • 10^'^ Mq is at best an upper 
mass limit of NGC 4449. A better guess is the dynamical 
mass inside the central gaseous ellipse which seems to be 
unaffected by interaction. E.g. the dynamical mass inside 
11 kpc would be a factor 3 smaller than Bajaja's value. 

Even more difficult is the mass determination of DDO 



the rotation curve exceeding the range observed by TuUy 
et al. by a factor of two. Moreover, they found no hint 
of a flattening or turn-over in DDO 125's rotation curve. 
This gives a lower mass limit of 4 • 10^ M©. Combining the 
mass estimates yields a lower limit of the mass ratio of 6% 
when using Bajaja's mass determination of NGC 4449. A 
more realistic mass ratio might be derived from the ranges 
of rigid rotation resulting in g « 18%. Thus, observation- 
ally it cannot be excluded that the mass ratio is below 
the critical one for the interaction scenario, however, the 
more realistic estimates seem to be in agreement with the 
values found in the simulations. The situation could be 
clarified when a reliable mass determination of DDO 125 
is available. 



4.2. Origin of the extended Hl-distribution around 
NGC 4449 

Another important question is related to the initial con- 
ditions of the previous simulations, i.e. the origin of the 
extended HI distribution. In principle, there are several 
possibilities: 

Idea 1. The HI comes from DDO 125. In a previous en- 
counter (e < 1) the gas was stripped off from DDO 125. 
Due to the long orbital period of e.g. 3.6 Gyr for e = 0.5 
the gas has enough time to settle down in a regular disc- 
like structure. However, the amount of HI gas in DDO 125 
is much smaller than the gas seen around NGC 4449. It 
seems unlikely that DDO 125 could have survived the loss 
of such a large fraction of its gas content without disrup- 
tion. 

Idea 2. The HI comes from NGC 4449. Though a star 
forming dwarf galaxy drives a galactic wind e.g. by means 
of type II supernovae, a detailed fine-tuning of the energy 
injection would be required to prevent a blow-out of the 
HI gas, especially if it acquired enough energy to form 
such an extended system. Moreover, the counter-rotation 
of the HI outside and inside the optical radius of NGC 
4449 makes this scenario not very appealing. 

Idea 3. The HI gas stems from a previous minor merger 
with NGC 4449- Though a satellite galaxy can be dis- 
rupted by the tidal field of NGC 4449, it is very unlikely 
that it is already disrupted far outside the galactic center 
where the HI gas is detected. 

Idea 4- The HI originates from the collapse of a dwarf 
companion. A dwarf galaxy forms in the vicinity of NGC 
4449 similar to the formation scenario of compact ellipti- 



125. Observations by TuUy et al. ( |1978| ) gave a dynamical 
mass of 5 • 10® M© insid e the optical radius of DDO 125. 
However, Ebneter et al. ( 1987 ) reported a linear increase in 



cals suggested by Burkert (1994). By this it undergoes a 
strong collapse which leads to mass ejection due to violent 
relaxation. In principle, a mass loss of up to 40% of the 
companion could occur, providing enough gas to explain 
the HI. However, where is the dwarf today? Excluding an 
unlikely radial orbit, the dwarf's orbit could in principle 
decay by dynamical friction. However, due to the large 
distance to NGC 4449 the decay time should be at least 
several orbital periods, which exceeds the Hubble time. 
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Idea 5. Another major encounter. Although NGC 4449 is 
a member of the loose Canes Venatici group, there seems 
to be no close galaxy which has a radial velocity in agree- 
ment with the rotational sense of the gaseous disc in NGC 
4449. The only candidate is the ringed Sab- type galaxy 
NGC 4736, located about 5.5 degrees south-east of NGC 
4449. The projected distance of 350 kpc and the differen- 
tial radial velocity of Av — 101 km s~^ (NED) is compat- 
ible with a scenario in which both galaxies experienced a 
close encounter about 3.5 • 10^ years ago on a hyperbolic 



orbit (Kohle 119991). Placing NGC 4736 at a distance of 4.4 
Mpc, its total HI m ass is 4.2 ■ lO^Mg (Bosma et al. |l977 ; 
Mulder & van Dricl 1993 ), which is quite low for Sab- type 
galaxies (Solanes et al. |l996| ). Might it be that NGC 4449 
formed its halo out of gas stripped from another galaxy? 

5. Summary and conclusions 

Several N-body methods are combined in order to develop 
a method for the determination of the parameters of in- 
teracting galaxies. This method is applied to the HI dis- 
tribution of NGC 4449. In a first step, the fast restricted 
N-body method is used to confine a region in parameter 
space which reproduces the main observed features. In a 
second step, a genetic algorithm (also using restricted N- 
body calculations) is employed which allows, in principle 
for both, an automatic fit of observational data even in 
a high-dimensional parameter space and/or a uniqueness 
test of a favoured parameter combination (only the latter 
has been done in this paper). For a genetic algorithm one 
typically has to follow a population of (at least) 100 mem- 
bers for 100 generations in order to get a good fit, provided 
the data are sufficiently accurate. Missing single pixels do 
not inhibit the parameter determination, as long as the 
key features are included in the data. Since a typical re- 
stricted N-body simulation takes a few CPU-seconds on a 
Sparc workstation, the whole fitting procedure is finished 
after 3-6 CPU-hours. 

In the third step, the results of the previous steps 
are compared with detailed self-consistent N-body sim- 
ulations. In the case of NGC 4449, they show that the 
restricted N-body calculations are reliable models for this 
encounter. The comparison with the sticky particle mod- 
els demonstrates that the HI gas can be modeled without 
any restriction by a purely stellar dynamical approach, 
provided the encounter is weak, as in the case of NGC 
4449 and DDO 125. 

From the previous simulations, we conclude that the 
extended HI features observed in NGC 4449 are created 
by an encounter with DDO 125. Prior to the encounter, 
the HI gas formed an extended, almost homogenous disc 
with a large scale length of about 30 kpc and a radial ex- 
tension of 40 (±10) kpc. The orientation of the disc is in 
agreement with the orientation of the inner ellipsoidal HI 
distribution. The orbital plane has an inclination angle of 
about 40° (±10°). The eccentricity is in the vicinity of 
a parabolic encounter, but favouring the region of ellipti- 
cal orbits (0.5 < e < 1.5). The apocenter distance of the 



galaxies is about 25 (±5) kpc. The closest approach hap- 
pened 3.5 — 6.2 • 10^ yr ago (depending on eccentricity). 
The mass ratio of both galaxies must be approximately 
0.2 (and cannot be smaller than 0.1). 
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Appendix A: The genetic algorithm 

In the following sections we briefly describe the main in- 
gredients of our genetic algorithm. The GA used here is a 
slightly modified version of P. Charbonneau's code pikaia 
(for details see Charbonneau ( 1995 )). The interface be- 
tween the GA and the simulations, i.e. the fitness determi- 
nations, is performed in a similar way as in Wahde ( 1998| ), 
however the coding of genes as well as the choice of pa- 
rameters used as genes differ from Wahde's approach. 



INTENSirf: original 



best fit ot generation 0007 



MAXiivlUM FiTNESS 




best fit of generation 0007 



FiTNESS DiSTRiBUTiON 





-20 20 



-20 20 



0.03 0.035 0.04 0.045 



Fig. A.l. Fitness and best fit of an early generation during 
the GA fitting procedure. The figure displays the artificial 
original data (left column: positions of the particles in the 
x-y-plane (upper left) and the resulting intensity contour 
lines (linear spacing, lower left), the best fit of the actual 
generation (middle column) and information about the fit- 
ness history, i.e. the maximum fitness (upper right) and 
the actual fitness distribution (lower right). The underly- 
ing grid in the left and middle column is used to determine 
the intensities. 

A.l. Fitness and parent selection 

In order to determine the fitness of an individual (i.e. a 
special parameter set), the N-body simulation has to be 
performed and compared with the observation. Generally, 
this is done by mapping the particles on a cartesian grid. 
However, here a quadratic grid of size 60 kpc centered on 
NGC 4449 is used. The resolution of the grid was set to 
7x7 in order to allow for a better resolution of the tidal 
features of an encounter. For each grid cell the intensity 
was calculated by summing the masses of the individual 
particles in that cell. The conversion to observed intensi- 
ties can be performed by assuming a mass-to-light ratio 
for the particles (or treating it as another free parameter) 
or by normalizing the intensities. Here the intensities are 
normalized to the total intensity of the whole grid. The 
quality can be measured by the relative deviation S of the 
intensities in both maps: 



E 

i{cclls) 



max(/rcf,j, /mod,i) 



/ 



The fitness / is then estimated by 
1 



1 



(A.l) 



(A.2) 



In the case of an impossible parameter configuration (e.g. 
a circular orbit with a distance less than the observed pro- 
jected distance) or any convergence problems the fitness 
of that point in parameter space is set to zero, prohibiting 
any offspring. Typically, only a few models in 10^ param- 
eter sets fail, most of them in the first few generations 
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before the system starts to approach the parameter space 
region of interest. 

In order to determine the parents of the next gener- 
ation, the members of the actual generation are ranked 
by their fitness. The parents are then selected randomly, 
whereas the probability of being chosen is proportional to 
the rank (roulette-wheel selection). 

A typical example of the st atus of the GA after a 
few generations is shown in Fig. |A.l| . In the left column 
the positions of the original data are discretized on the 
grid, yielding the intensities displayed on the lower left. 
Though the resolution of the grid (about 8.6 kpc) exceeds 
the thickness of the tidal structures, the large scale pat- 
tern is quickly found: in the middle column the best fit 
of the actual (here 7th) generation already agrees qual- 
itatively with the original. This demonstrates that large 
scale features already present in intermediate resolution 
maps can be sufficient to constrain the dynamics. This is 
in agreement with Wahde's (1998) calculations for equal 
mass galaxies: he recovered orbital parameters by using a 
4x4 grid. The quality of our fits is calculated only by a 
comparison of the lower two intensity diagrams. The fit- 
ness distribution in the population resembles a Gaussian 
at the beginning of the fitting procedure. 

A. 2. Coding of parameters 
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Fig. A. 2. Schematic diagram for the coding of parameters 
on a numerical chromosome. 



In organic systems the information is coded in genes in 
the form of sequences of the four nucleotide bases where 
a triplet of them (a codon) encodes one amino acid, the 
building block of proteins. The genes themselves are ar- 
ranged in one or several chromosomes which contain all 
the genetic information (or the genotype). The phenotype 
which is subject to the selection process is the result of 
the genetic information and the interaction with its envi- 
ronment. 

Translated to our fitting problem, the phenotype is the 
final intensity map created by a special choice of param- 
eters or initial conditions. The latter correspond to the 
proteins which physically express the genetic information. 
The mapping between the parameters and the genes or the 
number of chromosomes is not fixed in GAs. Many GAs 



apply a binary alphab et to encode the parameters (e.g. 
Holland 1975 ; Goldberg 1989). In Charbonneau's program 
PIKAIA a decimal encoding is implemented. For our calcu- 
lations we combined four decimal digits giving a number 
between zero and one. This gene has been mapped to a real 
parameter by a linear or logarithmic transformation, e.g. 
the orbital inclination is mapped linearily to the range [0° , 
180°] (Fig. [A.2| ). In general, the linear mapping of param- 
eter X, which is allowed to vary between a;iow and Xhigh, is 
performed by 



X — X\ovj -\- g ■ (^^high -^low) 



(A.3) 



Here g denotes the value of the gene which is in the range 
[0,1]. Similarly, the logarithmic mapping is done by replac- 
ing X by log(a;). All the genes are assumed to be located 
in a linear fashion on a single chromosome. Unlike organic 
evolution, the phenotype is completely determined by the 
genotype, no environmental influence on the individual 
acts after its genotype is fixed. 

A.3. Cross-Over 
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Fig. A.3. Schematic diagram for the cross-over operation. 



The main difference between evolutionary strategies 
and genetic algorithms is the sexual reproduction in the 
framework of GAs. This means that the chromosomes of 
the parents are combined in a new manner by means of a 
cross-over mechanism in addition to mutation. The prin- 
ciple idea is illustrated in Fig. A. 2 : After determining the 
two parents a random position on the chromosomes is se- 
lected. At this position the chromosomes are split into two 
fragments which are exchanged among the chromosomes 
resulting in two new ones. 



A.4. Mutation 



In order to introduce completely new information into 
a genepool a mutation process, i.e. a random change of 
the information on the chromosomes is introduced (Fig. 



A.4| ). For each position on a chromosome, mutation is ap- 
plied with a small probability (typ. Pmut ~ 0.5%). After 
selecting a position for mutation, an integer random num- 
ber in the range [0,9] replaces the original value at that 
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Fig. A. 4. Schematic diagram for the mutation process. 



position. Thus, a real change in the genetic information 
takes place with a probability of 90%, even if that position 
of a gene is subject to mutation. The overall probability 
for at least one real mutation on a chromosome coding 
6 parameters by 4 digits is about 11%. Thus, cross-over 
is the most important process introducing diversity into 
the population, whereas mutation affects only each 9th 
member of the population. In principle, the mutation rate 
can be increased, but the price would be a more or less 
random search in parameter space which is already after 
a few generations inferior to the oriented search of a GA 
method ( Wahde [T998| ) . 

A principal danger of GAs is inbreeding, i.e. the ho- 
mogenization of the genepool by one or a few dominant 
genes. If they successfully replace other genes, cross-over 
does not introduce sufficient (or - if inbreeding is complete 
- any) diversity into the population and thus any further 
improvement of the chromosomes is prohibited. The only 
way to overcome this situation is to introduce more diver- 
sity by mutation. Therefore, PIKAIA allows for a variable 
mutation rate which increases the mutation probability 
Pmut whenever the fitness distribution becomes narrow, 
i.e. inbreeding is suspected to operate. On the other hand, 
Pmut is decreased if the fitness distribution becomes very 
broad, i.e. the randomization is dominant. 



A. 5. Determination of the next generation 

After A^pop children have been created, the old generation 
is completely replaced by the next generation (full gener- 
ational replacement). The only exception is the survival 
of the fittest member of the parent generation, if it is fit- 
ter than the child which would replace it (elitism). This 
mechanism prevents the system from forgetting successful 
solutions. 

The initial population is randomly chosen in the ac- 
cessible parameter space. A standard population size is 
about Npop = 100 — 200 members. The GA converges 
typically after about TVgen = 100 generations. Thus, 10^ 
or a few 10** N-body simulations are required for a sin- 
gle GA fit. For a genetic algorithm simulation with 900 
particles (and no offset for an inner edge) 2 CPU-seconds 
are required per simulation compared to 3 CPU-hours on 
a GRAPES. Thus, self-consistent N-body models, such as 
the mentioned GRAPE models, would still need 3.4 years, 
whereas the restricted N-body models reduce the CPU- 
requirement to 5.6 hours on a ISOMHz-SparclO. 
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